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We derive several upper bounds for the superfluid stiffness D s for Bose and Fermi 
systems in terms of expectation values of local operators using linear response the- 
ory and variational methods. These give insight into the non-trivial dependence of 
D s on parameters such as disorder and interaction in systems with broken contin- 
uous translational invariance. Our best variational bound for disordered systems is 
obtained by allowing the phase twist applied at the boundary to be distributed in- 
homogeneously within the system. Path integral quantum Monte Carlo simulations 
are used to quantitatively compare the bounds and D a for disordered interacting 
Bose systems. 

I. Introduction: 

The superfluid density n s , or the closely related stiffness D s (defined below), characterizes the 
phase rigidity of a superconductor or a superfluid. Experimentally n s is more easily measured 
than the order parameter since external probes, such as magnetic fields or rotation, can couple to 
the relevant correlation function. Further, even from a theoretical point of view, n s is arguably 
the more fundamental characterization: e.g., in two dimensional systems n s can be non-zero even 
though the order parameter vanishes at finite temperatures. 

In translationally invariant (one component) superfluid systems, it can be shown on general 
grounds that at T = 0, n s = n, the total density: all of the fluid participates in superfiow [Q. This 
follows from translational invariance and the absence of low lying excitations that can couple to 
a transverse probe. However, as soon as one destroys continuous translational invariance, cither 
by the presence of a lattice or impurities, n s is no longer constrained to be the total density n, 
and becomes smaller as we shall describe below. The dependence of n s at T = on disorder is of 
considerable interest since with increasing disorder one could suppress the superfluid density and 
eventually drive the system into a non-superfluid, possibly insulating state 

In this paper we will derive exact upper bounds on n s (or D s ) which help understand how the 
superfluid density varies with interaction and disorder in both Fermi and Bose systems. An outline 
of the paper and a summary of our main results is as follows. Section II introduces the notation 
and definitions used in the rest of the paper. In Section III we discuss bounds for D s obtained 
from linear response theory, emphasizing the effect of breaking continuous translational invariance. 
Two simple upper bounds, eqn. ( p~0| ) and eqn. (|l^) , are derived and compared with results of path 
integral quantum Monte Carlo (QMC) simulations for non-disordered lattice boson problems. 

In Section IV we derive variational upper bounds on D s for a system with disorder by construct- 
ing a trial wave function for a system with twisted boundary conditions given the exact ground 



state for the untwisted case. We obtain two bounds: one, D° given in eqn. (21 ), which is the lattice 
analog of a result derived earlier by Leggett Q] for continuum systems, and an improved variational 
bound D* given in eqn. (Uq). The calculation of D* s is shown to be equivalent to calculating the 



conductance of a resistor network problem as explained in the discussion following eqn. (18). We 
also briefly discuss the generalization of these results to finite temperature. 
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The physical idea underlying the bounds D° and D* is the following. In a homogeneous system 
the applied twist at the boundaries is distributed uniformly across the entire system. However, 
in an inhomogeneous system the twist is accommodated preferentially in those regions of space 
where it costs the least energy, leading to a lower stiffness. Our results can be summarized by the 
following set of inequalities 

D a < D* < D° < tt(A'i) < Itixn (1) 

where D s is the exact supcrfluid stiffness, D° and D* are the bounds mentioned above, (K\) is 
(the absolute value of) the mean kinetic energy (KE) on a bond, t is the hopping strength in the 
lattice model and n is the filling or the average particle density. 

In Section V we illustrate the quantitative usefulness of our results by discussing various examples 
and comparing the bounds with D s obtained from QMC calculations on disordered Bose systems. 
(The details of the QMC are given in the Appendix.) The bounds are found to capture the non- 
trivial parameter dependence of D s , including its nonmonotonic variation in some cases. However, 
these bounds, which are based on local correlation functions, are not in general sensitive to the 
superfluid-insulator transition. In Section VI, we make several remarks about our results. Amongst 
other things we briefly discuss charge vs supcrfluid stiffness, describe a model for which D* captures 
a percolation transition that is missed by the other bounds, and obtain bounds on T c for 2D systems. 
We summarize our results in Section VII. 



II Definitions and Notation 

We consider interacting Bose or Fermi systems with a generic Hamiltonian 

r,a — l...d 

We work on a hypercubic lattice in d dimensions with lattice spacing a = 1, r runs over all N = L d 
lattice sites and a denotes the unit vector along the a th direction, c and are annihilation and 
creation operators; for Fermi systems we may also need spin labels which are omitted for simplicity. 
V includes all possible interactions between the particles as well as the possibly random (one-body) 
potential and is an operator which is diagonal in the number basis. Finally, t a (r) is the hopping 
amplitude between nearest neighbour sites r and r + a. This may again be a random function but 
is restricted to be positive. (Note that we do not consider frustration introduced by having both 
positive and negative hopping amplitudes.) 

The kinetic energy operator on the link (r, r + a) is defined as 

K a (v) = t a (r) [ct(r)c(r + a) + h.c] . (3) 

For later convenience, this definition differs from the standard one by a negative sign. The para- 
magnetic current along a at point r is given by 

j a (r) = -it a (r) [J(r)c{r + a) - h.c] . (4) 

The (paramagnetic) current-current correlation function for an iV — L d site system at zero tem- 
perature and zero frequency is given by 



where \tp n ) is the n th eigenstate of the system with energy ui n (n = denotes the ground state) 
and q is the momentum. The overbar denotes an average over the quenched disorder. (This 



is a shorthand notation for first evaluating the correlation function in real space, then disorder 
averaging and finally Fourier transforming to momentum space.) Now, within the Kubo formalism, 
we can define the superfluid stiffness D s in terms of the transverse current-current correlation 



— = {ih\Ki\ipo) ~ Xi (ft = 0,qx 0,u; = 0) (6) 

7T 

where = {<72, ■ ■ ■ , <Zd} is the set of momentum coordinates orthogonal to the 1 direction. The 
first term (Ki) = L^ d ^2 r ('tpo\Ki(r)\ipo) comes from the diamagnetic part of the response 0. 

An equivalent definition of the superfluid stiffness is arrived at by using twisted boundary condi- 
tions (b.c.) on the ground state wave function |^,^]. This can be imposed in the standard way by 
considering the system as a "ring" in the 1 direction, threaded by a flux (Note that $ = corre- 
sponds to periodic b.c.'s). In the site-occupation basis, moving one particle around the ring, keeping 
the others fixed, amounts to changing the configuration from {n} = {n(ri), . . . , n(r m ), . . . , n(r^)} 
with n(r m ) > 1, to {n'} = {n(ri), . . . , n(r m ) — 1, ... , n(rjy); n(r m + LI) = 1}. The twisted b.c. is 
then given by 

<{n'}|V*>=e**({n}|^). (7) 
The increase in ground state energy E((fr) = (i]j^\H\ip^) for a small applied twist 

E{*)-E$) = \^L d (j^ + ... (8) 

defines the superfluid stiffness D s . This is related to the superfluid density n s and the helicity 
modulus T S via 

5£ = £ «,i m ^)= 2t n. = r. (9) 

Note that for a disordered system, the twist definition of D s is valid for any given disorder real- 
ization, and thus disorder averaging can be done at the end of the calculation. 

The generalization of the above definitions of D s to finite temperatures is straightforward. In 
the Kubo formula (||) and (^) , the ground state expectation value has to be replaced by a thermal 
average, and correspondingly in the twisted b.c. formulation (||) the energy has to be replaced by 
the free energy. 

Ill Bounds from Linear Response Theory 

From (§ we see that xi(q, w = 0) is non-negative. It then follows from that 

D s /n < (Ki). (10) 

At finite temperatures the right hand side should be viewed as thermal expectation value. 

This bound can be further simplified for the case of uniform hopping t a (r) = t. Using the 
Cauchy-Schwarz inequality and the fact that the geometric mean is less than the arithmetic mean 
we find 

(iCi(r)) = <((c f (r)c(r + 1) + h.c.)) < 2tyjn(r)n(r + 1) < t(n(r) + n(r + i)), (11) 
where n(r) = (c^(r)c(r)). Further (K\) < tL~ d J2 T ( n ( r ) + n ( r + 1)) = an d using (0) we get 



n s < (K 1 )/2t < n. 



(12) 



Thus, the superfluid density for an interacting lattice or disordered system can be, and in general 
is, less than the average particle density. The aim of this paper is to understand how this happens 
and (in the following Sections) to derive bounds that improve upon the ones above. In the present 
section we focus on the non-disordered problem; we turn to the the disordered problem in later 
sections. 

Let us begin by emphasizing the importance of broken continuous translational invariance for 
n s (T = 0). For a translationally invariant system, the diamagnetic piece (the analog of (-Kj.) in 
(((])) is given by n. Furthermore the q = current operator is the total momentum which commutes 
with the Hamiltonian, leading to a vanishing numerator in (|^). Thus for a continuum superfluid 
Xi vanishes and n s (T = 0) = n since there are no low lying excitations which can couple to a 
transverse probe. (The only low lying excitations present are longitudinal phonons, and even these 
may be gapped for a charged system). However in the presence of a lattice or impurities, the total 
(q = 0) current operator does not commute with the Hamiltonian in general, thus xi 1S non-zero, 
and ( |10|) is in fact a strict upper bound on the stiffness. 

It is important to note that approximate treatments can miss out on the paramagnetic current 
correlation piece and give just the diamagnetic part even on the lattice. For example, the BCS 
mean field theory for superconductors reduces the problem to one of non-interacting particles 
(Bogoliubov quasiparticles) with a gapped spectrum and a current operator which commutes with 
the reduced BCS Hamiltonian (even on a lattice). Thus, estimating the reduction of the superfluid 
stiffness from its diamagnetic value at T — requires including fluctuations beyond mean field 
theory @. 

We next discuss interacting (non-disordered) Bose systems on a lattice to illustrate the usefulness 
of the bound (|l0|). We consider the 2D boson Hubbard model with Hamiltonian 

H hosc = -t J2 W(r)c(r + a) + h.c.] + ^ ^ n(r) [n(r) - 1] (13) 

r,Q = l. ..d r 

where n(r) = c^(r)c(r). By comparing the bound with exact results from QMC simulations on 
finite two-dimensional lattices, we gain insight into how the superfluid stiffness for lattice systems 
depends on the interaction, unlike in the continuum case. The calculation of D s from QMC is 
described in the Appendix. Here we summarize relevant numerical details. The QMC simulations 
have been done on lattices of size N = 6 x 6 at commensurate and incommensurate fillings with 
filling fractions n = 1 and n = 1/2 respectively. The inverse temperature (3t = 4 is discretised into 
N T — 32 steps which implies that St = (3 /N T — 1/8 is the elementary time step on which the Trotter 
approximation is done. We averaged over 10 s MC sweeps through the space-time lattice with 10 6 
equilibration sweeps during which we did not collect data. The very large number of sweeps was 
necessary to get small enough error bars for the superfluid stiffness D s in order to distinguish the 
QMC results from the bounds, especially for large interaction and (or) disorder which is the region 
of parameter space explored in both this as well as later sections. The calculations for each data 
point took approximately 36 hours on a 100 MHz pentium. 

In all figures, here and below, we will plot the superfluid stiffness normalized such that it is unity 
for the noninteracting clean (no disorder) problem at T — 0, independent of the filling factor. We 
normalize the bounds also by the same factor. Thus, at T = 0, what we plot for the stiffness is 
really n s /n as can be seen from the definitions above since n s = n for the clean noninteracting 
problem at zero temperature. 

First consider the case of incommensurate (non-integral) filling where the system is expected to 
be superfluid for all values of the interaction strength < U < oo. We see from Fig. 1 that the 
superfluid stiffness D s , though non-zero for all U, decreases monotonically with the repulsion U. 
The bound n(Ki) derived in eqn. (|l0|) tracks the behavior of D s fairly well even at large U . The 



error bars on the bounds are much smaller than those on D s and are hence not indicated in Fig. 1 
or any of the subsequent figures. 

Next, consider a commensurate filling of one boson per site (on average). With increasing U we 
expect a Mott transition Q from a superfluid into a Mott insulator at a critical value U = U c . 
Our aim here is not to locate and characterize this transition very precisely (which would require 
a finite size scaling analysis of the QMC data as in Ref. pi]]). In Fig. 2, we compare the KE bound 
with D s obtained from QMC on the finite lattice. Even though the bound captures the overall 
decreasing trend of D s , it falls off as t 2 /U for large U (as expected from perturbation theory) 
rather than vanishing like D s for roughly U > 20. Being a purely local quantity the KE bound 
misses the transition. This brings out the importance of the excitations which contribute to Xi 
and reduce the stiffness from its diamagnetic value eventually driving D s to zero at the transition. 

It hardly needs to be emphasized that any non-zero upper bound on D s such as (fio|), though 
valid, will not be useful for systems which are in fact non-superfluid (metallic or insulating). 
In such cases the paramagnetic piece xi ignored in the bound ( |l0| ) must exactly cancel (-Ki), 
leading to D s = 0. The most trivial example is the non-interacting fermi gas: the total current 
operator commutes with the Hamiltonian even on a lattice, leading to a vanishing numerator in 
the paramagnetic term in (^) , but the low lying particle-hole excitations contribute in such a way 
that there is complete cancellation between the diamagnetic and paramagnetic terms in (j^) . More 
non-trivial examples are the Mott phase of the boson Hubbard model discussed above, as well as 
the large U (fermion) Hubbard model at half-filling, which is also a Mott insulator. In both cases, 
the total current operator no longer commutes with the Hamiltonian and the high energy states 
(excitations across the charge gap) contribute to xi leading once again to a complete cancellation 



IV. Variational Bounds 

We now turn to the task of improving upon the bounds of the previous Section, especially for 
disordered systems described by the Hamiltonian in eqn.(^) for an interacting Bose or Fermi system 
in the presence of a site-dependent disorder potential. To this end we use a variational method 
first suggested by Leggett jl); see also |fL2| . Although these authors worked in a first quantized 
representation, we find it more convenient to work in the occupation number basis. Let us assume 
that the exact (normalized) ground state \4>o) of the Hamiltonian (^) with periodic boundary 
conditions (b.c.) in all directions is known. We make an ansatz for the ground state of (|j) 
subject to the twisted b.c. (Q) along the i direction, with periodic boundary conditions in the 
other (d — 1) directions. Let: 



|t/>*} = exp 



;^0(r)n(r) 



l^o) (14) 



where r runs over all lattice sites and 9(r) are the variational parameters at our disposal. The 
twisted b.c. (Q) imposes the constraint 

0(ri=L+l i T±)-9(r 1 =l,r±) = §. (15) 

where is the coordinate in the 1 direction along which the twist is applied and = r 2 , . . . , 
refers to the set of coordinates transverse to this direction. The variational estimate for the 
difference in energy between the twisted and the untwisted boundary condition cases is then (for 
small twist); 

= \ E (0(r)-e(r + a)f(K a (r)) (16) 

r,et=l...c£ 



where (K a (r)) — (ipQ\K a (r)\ipo) . ( ln this expression the term linear in [0(r) — 9(r + a)] vanishes, 
since we assume that the ground state |^o) does not break time-reversal symmetry and is thus real 
and carries no current on any link.) Minimizing Atrial with respect to the parameters 9(r) (for 
those r not on the boundary on which the twist is applied) we obtain the following set of equations: 

]T {(K a (r)) [0(r) 0(r + d)] + <#_ a (r)) [0(r) - 9(r a)}} = (17) 

a=l...d 

where we have defined K- a (r) = i_ Q (r) [c^(r)c(r — a) + h.c] in parallel with (||). 

There are in all L d — L d ~ x equations implied by eqn.(|T7|) and L d_1 equations implied by the 
constraint eqn.(|l5|) thus giving us a total of L d equations which allows us to solve for the phase 
9*{y) at each point of the lattice. The solution {9*(r)} of equations (J!?]) and ( ]l^ ) defines the 
optimal choice of the local phases in in response to the applied twist. This in turn determines 
.D* via the relation 

l^L d - 2 ($) 2 = Atrial = \ E (O*(r)-0*(r + a)) 2 (K a (r)). (18) 

r,a— l...d 

It follows from the variational principle that AE* lial [<S>] > AE[$] and thus D* > D s from eqn. (|l8|) 
above. 

It is not possible to start from the set of equations ([l7]) and (|l5|) , and arrive at an analytic closed 
form solution for the bound D* in general in d > 1. (Later we will discuss the d = 1 case and also 
a simplified form which admits an analytic solution in arbitrary dimension) . It is however possible 
to calculate D* numerically solving the linear equations ( filf ) and (|l7|). For this, we first map the 
problem of obtaining D* on to a particular random resistor network problem which gives us some 
intuition for the results. The lattice points r of the above system form the nodes of the network, 
(r) is the voltage at node r, and {K± a (r)) is the conductance of the link (r, rid). The equations 



(17) represent (Kirchoff) current conservation at each node r, and the constraint (|15|) is a constant 
voltage $ maintained across this resistor network. The quantity to be minimised in eqn. ( ^6|) is 
the (Joule) heat dissipation in the network. Thus L d ~ 2 D*/n, defined by (|l^), is the conductance 
of an equivalent resistor which dissipates the same amount of heat as the network for the same 
voltage difference across its ends. 

It is instructive to write down the solution to the above equations in d = 1 where a closed form 
result can be obtained. We find 

(^))E»=i...l (Ki(x)y 



>( x + l)-9*( x ) = 7 — 7 -^- —————^ (19) 



from which we get D* = LitK^ot with 

^ = S..L^^y (20) 

It is easy to see that 1/JCtot is simply the total resistance of the 1— d network obtained by adding the 
resistances of the links in series. In order to help visualise this result, we show exact diagonalization 
results for a system of bosons on a finite one dimensional lattice in Fig. 3. The Hamiltonian used 
is Hbose + X) r ^( r )'^( r ) where Hb ose is defined in cqn.(|l3[). The parameter values chosen are given 
in the figure. We plot the way the imposed phase twist at the boundary is accommodated along 
the chain, and find that a large phase twist occurs across bonds where the local kinetic energy is 
small as expected from eqn. (pj|). 

The d = 1 result suggests a simple generalization to arbitrary d. If we assume that 0(r) depends 
only on r%, the coordinate along the applied twist direction, and is independent of the transverse 



coordinates rj_ = {^2, --j^d} then (following the steps in the derivation of (filf)) we get an explicit 
analytic expression for the bound D° as 

1J » < u * — r2-d i Con 

t ~ t EntErJ^iCri.rx))]- 1 1 ' 

where (K±(ri, rj_)) = (^(r)). We note that D° is the lattice analog of the continuum result 
obtained by Leggett Q and Chester [^2[ ; in the continuum Kx(r) gets replaced by the density 
n(r). The two bounds D* and D° coincide in one dimension as seen from ( po|) and (^l|). The 
two bounds also coincide for the non-disordered problem in any d; the (discrete) translational 
invariance ensures that the bound is just the mean kinetic energy on the links in the direction of 
the applied twist. However, in general, for d > 1, D* < D° and is a better bound because it has 



more variational parameters |13|. 

We have thus derived all the inequalities in ([!]), except for the third: D° < k{K\). The first two 
inequalities were derived above and the last one in ([!]) , valid only for systems with uniform hopping, 
was derived in ([l2|). To complete our task, we use the inequality M~ x J2i a i ^ [M _1 J3j a^ -1 ] -1 
for positive numbers a;, i = 1, . . . , M, in (pTf). This yields 

? ^"fEij^OT'™- <22) 

Generalization to Finite Temperatures 

The above bounds can be easily generalised to finite temperatures using a trial (variational) 
density matrix as sketched below. This can be easily shown to lead to the same equations as 
the zero temperature case with the only difference being that we get thermal expectation values 
instead of ground state expectation values for the local kinetic energies. This allows us to compare 
the bounds with the QMC results for D s in the next section. The local kinetic energies needed for 
calculating the bounds are also obtained using QMC as described in the Appendix. 

Let po — e~P Em \m)(m\ be the density matrix for the untwisted case, where |m) denotes 
the eigenstate of energy E m for the system with untwisted boundary conditions. Consider a trial 
density matrix for the twisted case of the form 

/5($) = exp[i^n(r)0(r)] p exp[-i^n(r)0(r)]. (23) 



with boundary conditions as in eqn(|15[). The trial free energy is given by 

F(Q) < F W (*) =Tr[p($)£] + p- 1 Tr [p($)lnp($)] (24) 
and we need to minimize -Final (*&) as before. Using the unitary nature of the transformation in 



eqn.(23) and the cyclic property of the trace, it is straightforward to show that the entropy term 
/3 _1 Tr[p(<i>)lnp(<i>)] remains unchanged even when the twist is put in. The energy term Tr (p(<f>)H) 
leads to an expression identical to that obtained at T = with thermal expectation values instead. 

V. Comparison with Quantum Monte Carlo Results: 

In this Section we illustrate the usefulness of the bounds by comparing them with the stiffness D s 
obtained from QMC [ p~4|JT~^ ] for interacting Bose systems described by the Hamiltonian in eqn.(||) 
with t a (r) = t and 



F= ¥ ^ri(r)(n(r)-l) + ^nr)ri(r) 

r r 



(25) 



We work in two dimensions and choose the site dependent (disorder) potential V(r) as given below 
for various cases. 

We first consider a particular case of broken (discrete) translational invariance to illustrate the 
difference between the various bounds in ([!]) . We choose a certain set of sites on the lattice where 
V(r) — V > (marked in the inset to Fig. 4) while fixing the potential at other sites to be zero. We 
then study the effect of varying V on the superfluid stiffness D s and compare it with the various 
bounds, all computed for a twist applied along the x— direction. The results are shown in Fig. 4: 
as V is increased, the mean KE bound and D° are both considerably larger than the true stiffness 
while the bound D* tracks the true stiffness extremely well. For instance at large V , as seen from 
Fig.4, (Ki) ~ 0.8, D°/2wnt ~ 0.5, while D*/2irnt and D s /2irnt are both - 0.05 and thus D* 
is about an order of magnitude better than the other bounds. We do not know if the excellent 
agreement between D* and D s would persist in the thermodynamic limit, nevertheless it is clear 
that D* does considerably better than the other bounds. 



We next consider the disordered Boson Hubbard model, which is defined by (25), with V(v) 
at each r an independent random variable uniformly distributed in the interval [— Vdis> Klis]- We 
study the case of incommensurate filling, i.e., non-integral average site occupancy. As mentioned 
earlier, the bounds are valid for any given disorder realization. We therefore choose to consider a 
particular realization of the random potential in Figs. 5 and 6 in order to illustrate the bounds. 

In Fig. 5 we keep the disorder strength Vdis fixed, and look at the interaction [/-dependence 
of the superfluid stiffness and the bounds. For small U, we see that D s is small, since at U = Q 
we expect that the system is in a localized phase, with all bosons in the lowest lying, localized 
single-particle eigenstate. Note that due to the finite size of the calculation, D s does not vanish 
even at U — 0. Nevertheless, the bounds D* and D° do capture the large suppression of D s at 
small U, while the mean kinetic energy is clearly unable to do so. For intermediate U we find 
nonmonotonic [/-dependence in D s fI3j . From a mathematical point of view, there is no guarantee 
that because D s is a nonmonotonic function of some parameter, say U , an upper bound on D s 
should show similar [/-dependence. It is therefore interesting that this feature is captured by all 
the bounds. At larger values of U, D* and D° are closer to the KE in value, than to D s , but all 
three track the overall trend of D s fairly well. 

In Fig. 6 we keep the interaction U fixed and study D s and the bounds as a function of disorder 
strength Vdi s - We find that D s decreases monotonically with disorder. This trend is captured by 
all three bounds, however, we clearly see that D* does increasingly better than the other bounds 
in the large disorder regime. 

VI. Remarks: 

In this Section we collect together some remarks which are meant to clarify different points which 
we feel are important. We also consider a model where the bound D* captures a transition unlike 
in the examples presented earlier. 

Charge Stiffness versus Superfluid Stiffness 

A bound similar to the mean kinetic energy bound ( |lo|) has been discussed in Ref. [fl6f in connec- 
tion with the charge stiffness, which is the coefficient of 5(u>) in the optical conductivity of a perfect 
conductor. We emphasize that although the linear response formalism leads to essentially identical 
bounds for both the charge stiffness and the superfluid stiffness D s , the variational wavefunction 
construction of Section IV and the associated bounds are valid only for D s . Recall that the charge 
stiffness is defined in terms of the energy of the adiabatically generated ground state as the twist is 
introduced j(| . The variational ansatz that is considered here is by construction a trial ground state 
of the system at finite twist and is not an adiabatic continuation of the untwisted ground state. 
We therefore emphasize that unless there is a gap which makes the charge and superfluid stiffness 
equal fel, the above variational construction does not give us a bound on the charge stiffness. 



Kinetic Energy : Exact versus Approximate Estimates 

We should note that our bounds are expressed in terms of the local kinetic energies evaluated 
in the exact ground state which is in general not known. Nevertheless, one can derive useful 
qualitative information simply from the form of the bounds. In order to get quantitative estimates 
of the bounds and of D s one needs to resort to numerical methods like QMC. However, the bounds 
which involve local correlations, are much less affected by finite size errors. 

If a trial ground state wave function is used, instead of the exact one, to evaluate the local 
kinetic energies, then the expressions derived above do not lead to rigorous upper bounds. Such 
approximations could, nevertheless, lead to useful estimates. 

Quantum Percolation Model 

We now describe an example which highlights a qualitative difference between the bounds D* 
and D° and where D* in fact captures a phase transition unlike earlier examples of both clean and 
disordered systems. Consider a "quantum bond percolation model" of bosons on a lattice with 
hopping disorder such that the hopping matrix element between two neighbouring sites is either t, 
with probability (1 —p), or zero with probability p. We assume that the twist is applied along a 
certain direction and that we have fixed boundary conditions in the perpendicular directions. 

As the probability p of having a "broken link" increases, the superfluid stiffness reduces until 
it vanishes at some p = p\ oc via a (quantum) localization transition. We do not expect either 
D* or D° to show a signature of this transition. The bound D* would go to zero when there is 
any set of broken links that can break up the lattice into two disconnected parts across the twist 
direction. Thus we expect D* to vanish at p = p c corresponding to the classical bond percolation 
transition and p c > pi oc - In order for D° to vanish, the lattice needs to be disconnected along a 
(hypcr)plane perpendicular to the twist direction. Since the hopping matrix elements are chosen 
randomly, realizations with such a set of broken links will be of measure zero. Thus D° will not 
capture this percolation transition. 

Bounds on T c in 2D: 

In two dimensions an upper bound on the superfluid density allows us to put an upper bound on 
the superfluid to normal transition temperature T c in a clean system. For this we use the relation 
between the jump in superfluid density at the Kosterlitz-Thouless (KT) transition and the KT 
transition temperature @ D 3 (T~) = 2T C (with k B = 1). Thus, D S (T~) < D*(T C ) < D*(T = 0) 
which gives us the following bound on T c : 

T c < D* s (T c )/2 < D* S {T = 0)/2. (26) 

Note, however, that in general a lower superfluid density at zero temperature does not necessarily 
imply a lower T c . A counterexample is the case of s-wave BCS superconductors with weak disorder: 
Anderson's theorem tells us that in this case T c and other thermodynamic properties are unaffected 
but one would expect D s to be suppressed by the disorder. 

VII. Conclusions 

In this paper we have obtained several upper bounds on the superfluid stiffness at zero temper- 
ature as well as finite temperatures. These have been obtained using the Kubo linear response 
formalism as well as variational estimates. These bounds give qualitative insight into the depen- 
dence of D s on various parameters. 

For clean systems, the average kinetic energy (Ki) gives a good upper bound on D s /ir and pro- 
vides insight into the difference between lattice and translationally invariant (continuum) models. 
The bound captures the qualitative trend of the superfluid stiffness D s to decrease with increasing 
repulsive interaction U. For incommensurate filling, the bound compares well with D s at all values 



of the interaction. At commensurate filling, the bound is good at small values of U but fails to 
capture the superfluid-Mott insulator transition at large U. 

For disordered systems we obtain two improved variational bounds D* s and D° (D* < D° < 
■k{K\)) by allowing larger phase twists on links with lower local kinetic energy. We quantitatively 
compare the bounds with QMC results and find that they capture non-trivial parameter dependence 
of D s , including its nonmonotonic variation in some cases. 



Appendix 



Path Integral Quantum Monte Carlo 

In this Appendix we briefly describe the path integral QMC method used to calculate the 
local kinetic energies and the superfluid stiffness D s presented in the paper. For more details see 
Refs. [ p~8]JTg| ] . The notation in this appendix differs from that in the rest of the paper since here 
we work in first quantized representation. The position vector of the n th boson, n = 1 . . . Nj,, at 
imaginary time t, where < r < /3, is denoted by r„(r). The diagonal density matrix for bosons 
in the canonical ensemble is given by 

p(R, R;P) = ^ E^l e*p(-0H)\V[R]) (27) 

which describes the amplitude for the A^-particle system starting in a configuration R = 
{ri,r2, . . . ,rjv b } at time r = to be in configuration T^i?], where V is a permutation of R, at 
time t = [3 = 1/T. The partition function is given by Z = J2r p(R, R] /?)■ Each particle describes 
a path or a world-line in the r-r space, and the evaluation of the partition function amounts to 
summing over all possible paths which is most efficiently done by QMC techniques. 

The first step in evaluating the density matrix is to break the total time (3 into small steps 
St = j3/N T . Upon inserting complete sets of states we get 

p(R,R-,0) = ^J2zZ--' E p{R,Ri\St)p{R 1 ,R 2 -8t)---p{R (n ^ 1) ,V[R];5t) (28) 

V -f?i -R(jv T -i) 

A particular configuration of the particles on all the time slices, is denoted by 

S = [P; R, Ri,R 2 , ■ ■ ■ R(n t -i), Rn t ] (29) 

and looks like a collection of strings. Our aim is to sample the probability function 

Prob[5] = p(R, R 1 ;5t) P (R 1 ,R 2 ; St) ■ ■ ■ p(R^-i),V[R];St) (30) 

We use the Metropolis importance sampling method to generate a sequence of M string configu- 
rations Si , S%, . . . Sm such that in the limit M — > oo the probability of finding a particular config- 
uration 5^ in this sequence is proportional to Prob^^]. The actual implementation of the Monte 
Carlo method, including importance sampling, and the bisection method for efficient exploration 
of phase space, is described in Ref [pTofl . 

Calculation of the Local Kinetic Energy 

For simplicity, first consider an operator O = fjfij) which is diagonal in the position basis, 
like the interaction energy in the second term of eqn.(|Rj) (similar considerations apply to the 
disorder energy in eqn.(]25|)). Then 

^l— 1 r— 1 i 

where nf(r) is the number of bosons at site i on time slice r in the string configuration S^. 

Next, consider the local kinetic energy operator given in eqn.(||). Recall that in the Monte Carlo 
procedure, the density matrix in an elemental time St is approximated by 



A',, 



p(R T ,R T+1 ) = J] Pl «(r),r^(T + l)) x exp 



-(5r/2)J2[K(r)+Ut(r + l)} 



(32) 



where p\ is the single particle density matrix corresponding to the kinetic energy and disorder 
potential operators and ^(n^r)) = (f//2)n^(r)(n^(r) — 1) arises from the interaction between 
the bosons. 

To evaluate the kinetic energy, pick a string configuration with the n th boson at (r) on time 
slice r and imagine moving that boson to one of its four neighboring sites r^f (r) on the same time 
slice generating a new configuration. For clarity, we suppress the /i index labelling a configuration 
in the following equations. The change in potential energy is given by U new — U id- 

U old =U[n(r n ( T ))] +U[n{r' n ( T ))} (33) 

and 

U new = U Kr„(r)) - 1] + U Hv'Jt)) + 1] (34) 

There will also be a change in the single particle density matrix. 
Thus the local kinetic energy for a given string configuration is 

K(r n (r),r' n (r)) = ^Ml^L±llj exp[-£/„ etu + U M ] (35) 

which is then averaged over all the configurations selected by the Monte Carlo procedure. 

Calculation of Superfluid Stiffness D s 

The current along the a th direction is given by 

N„ 

ja(r) = ^[r„ Q (r + St) - r na (r)]/5r (36) 

n=l 

where r na (r) is the a th component of r„(r) (a = 1 . . . d). The winding number W a which counts 
the number of times the path of a particle winds around a box with periodic boundary conditions 
is given by (2^] 

1 Nt 

W a = w Y,Ur) (37) 

T=l 

and the superfluid stiffness is given by 



^-L^Y^iWl) (38) 



7T 



In the case of Fig. 4, we calculate the winding number, the stiffness and the bounds all along the 
x— direction only as indicated in the figure inset. For all other figures, we calculate the total D s 
given in eqn.(5|). 
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FIG. 1. The superfluid stiffness D s and the KE bound n{Ki) obtained from path integral QMC for the 
Boson Hubbard model in eqn. (13) at an incommensurate filling n — 0.5 for parameter values shown in the 
figure. The QMC data here (and in other figures) was averaged over 10 8 runs to sufficiently reduce error 
bars in order to allow comparison with the bounds. The normalization here and in the remaining figures 
is chosen to make the normalized stiffness unity for a clean non-interacting system at T = independent 
of the filling factor. 




FIG. 2. Superfluid stiffness D a and KE bound -k(Ki) for the Boson Hubbard model (13) plotted for 
commensurate filling n and parameter values indicated in the figure. In the thermodynamic limit, D s 
would vanish for U > U c (where U c /t ~ 20) signalling the superfluid to Mott insulator transition while on 
a finite system it remains non-zero even beyond that value and is zero (within error bars) only at large U. 
The bound however is clearly non-zero even at large U and can be shown from perturbation theory to go 
as t 2 /U for t >> U. The bound does not capture the Mott transition. 
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FIG. 3. The phase 6(r) calculated at each site (shown as arrows) indicating how the applied phase twist 
$ = 7r is accommodated along the lattice in a one-dimensional case for a particular realization of disorder 
for the dirty boson problem of eqn. (25). The absolute value of the kinetic energy is indicated along the 
links and the disorder potential is shown in the lower part of the figure (both in units of t.) The "bad 
bonds" (the first two links with smaller kinetic energy) allow for a greater phase twist across them as 
shown. The data was obtained by exact diagonalization for a system of 4 bosons on a 6-site chain. 




FIG. 4. Superfluid stiffness plotted for a 6 x 6 system at incommensurate filling in which the sites 
indicated in the inset are kept at a potential V with the rest kept at a potential of zero. The phase twist 
is applied along the direction shown by the arrow in the inset and the bounds and the stiffness D s are 
obtained only along that direction. The error bars on D s are of the size of the data point and hence not 
indicated. The difference between the various bounds is evident. Notice that the bound D* is better than 
the other bounds by about an order of magnitude. 




FIG. 5. Superfluid stiffness plotted as a function of interaction U for the dirty boson problem in 
eqn. (25) at incommensurate filling for Vdis/t = 10 and other parameter values indicated in the figure. The 
bounds capture the non-monotonic behavior of D s obtained from QMC as a function of U. The bound Dl 
is clearly better than the KE bound when compared with D s . At large U, the kinetic energies on various 
links are more uniform and the bounds are closer to each other. 




FIG. 6. Stiffness plotted as a function of disorder, Vdis, for the dirty boson problem in eqn. (25) at 
incommensurate filling for parameter values indicated in the figure. The bound D* gets progressively 
better than the other two bounds in the disorder dominated regime. This is precisely the regime where 
the local kinetic energies are very inhomogeneous. 



